"""Rate law factory and classes"""

import math
import numpy as N
import phyconst

class RateLaw():
    """Rate law base class"""
    def __init__(self, type, heatX, carbonNum):
        self.type = type
        self.mHrx = heatX[0]
        self.synHrx = heatX[1]
        self.carbonNum = carbonNum
    
    def calc_rate(self, pArray, temp):
        raise NotImplementedError('RateLaw.calc_rate')
    
    def calc_heat(self, rArray, sels):
        raise NotImplementedError('RateLaw.calc_heat')
        
class LH2RateLaw(RateLaw):
    """L-H-2"""
    def __init__(self, params):
        RateLaw.__init__(self, 'L-H-2', params[-2:], 1)
        [self.k10, self.E1, self.k20, self.E2] = params[0:4]
        
        [self.B1, self.C1, self.alpha1, self.beta1, self.gamma1,
         self.delta1, self.epsilon1, self.m1] = params[4:12]
         
        [self.B2, self.C2, self.alpha2, self.beta2, self.gamma2,
         self.delta2, self.epsilon2, self.m2] = params[12:20]
    
    def calc_rate(self, pArray, temp):
        """Calculate the rate law, where pArray is in Pa. It needs to be converted
        to bar
        """
        pCO = pArray[0] / 1e5
        pH2 = pArray[1] / 1e5
        if pCO < 1e-6 or pH2 < 1e-6:
            return N.array([0.0, 0.0, 0.0, 0.0])   
        rCO = self.k10 * math.exp(-self.E1 / temp / phyconst.Rgas) * \
            pCO ** self.alpha1 * pH2 ** self.beta1 / \
            (1 + self.B1 * pCO ** self.gamma1 * pH2 ** self.delta1) ** self.m1
        rCH4 = rCO * self.k20 * math.exp(-(self.E2 - self.E1) / temp / phyconst.Rgas) * \
            pCO ** self.alpha2 * pH2 ** self.beta2 / \
            (1 + self.B2 * pCO ** self.gamma2 * pH2 ** self.delta2) ** self.m2
        if rCH4 > rCO:
            rCH4 = rCO
        rH2 = 3 * rCH4 + 2 * (rCO - rCH4)
        rH2O = rCO
        return N.array([-rCO, -rH2, rCH4, rH2O])
    
    def calc_heat(self, ra, sels):
        rCO = ra[0]
        selCH4 = sels[0]
        deltaH = rCO * (self.mHrx * selCH4 + self.synHrx * (1 - selCH4))
        return deltaH
    
class DetailRateLaw(RateLaw):
    """Detail rate law in Yi-Nin Wang's paper"""
    def __init__(self, params):
        carbonNum = params[0]
        if carbonNum < 2 or carbonNum > 50:
            raise ValueError('Carbon number in Detail rate law needs to 1 < n <= 50')
        RateLaw.__init__(self, 'DETAIL', params[-2:], carbonNum)
        self.carbonNum = params[0]
        [self.k1, self.k5M0, self.E5M, self.k50, self.E5, self.k60, self.E6] = params[1:8]
        [self.kv0, self.Ev, self.k_6, self.Kv, self.K2, self.K3, self.K4] = params[8:15]
        self.temp = 273.15
        (self.k5M, self.k5, self.k6, self.kv, self.kp) = N.zeros(5)
    
    def calc_rate(self, pArray, temp):
        """pArray in Pa, it will be converted in bar"""
        pArray = pArray / 1e5
        pCO = pArray[0]
        pH2 = pArray[1]
        pH2O = pArray[3]
        pCO2 = pArray[4]
        pParaffin = pArray[5 : 5 + self.carbonNum - 1]
        pOlefin = pArray[5 + self.carbonNum - 1:]
        if pCO < 1e-6 or pH2 < 1e-6:
            return N.r_[[0.0, 0.0, 0.0, 0.0, 0.0], N.zeros(2 * (self.carbonNum - 1))]     
        if self.temp != temp:
            self.temp = temp
            self.k5M = self.k5M0 * math.exp(-self.E5M / temp / phyconst.Rgas)
            self.k5 = self.k50 * math.exp(-self.E5 / temp / phyconst.Rgas)
            self.k6 = self.k60 * math.exp(-self.E6 / temp / phyconst.Rgas)
            self.kv = self.kv0 * math.exp(-self.Ev / temp / phyconst.Rgas)
            self.kp = 5078.0045 / temp - 5.8972089 + 13.958689 * 1e-4 * temp - 27.592844 * 1e-8 * temp**2
        
        # alphaA: Anderson-Schulz-Flory(ASF) chain growth factor
        # betaN: olefin readsorption factors, N >= 2 carbon
        # alphaN: overall chain growth factors, N >= 1 carbon
        alphaA = self.k1 * pCO / (self.k1 * pCO + self.k5 * pH2 + self.k6)
        k11 = self.k1 * pCO / (self.k1 * pCO + self.k5 * pH2)
        k12 = self.k_6 / (self.k1 * pCO + self.k5 * pH2 + self.k6)
        a1 = alphaA ** N.arange(1.0, self.carbonNum)
        g = lambda n: sum(alphaA ** N.arange(0.0, n - 1) * pOlefin[0:n-1][::-1])
        a2 = N.array([g(n) for n in range(2, self.carbonNum + 1)])
        betaN = self.k_6 / self.k6 * pOlefin / (k11 * a1 + k12 * a2)
        alphaN = self.k1 * pCO / (self.k1 * pCO + self.k5 * pH2 + self.k6 * (1 - betaN))
        alpha1 = self.k1 * pCO / (self.k1 * pCO + self.k5M * pH2)
        alphaN = N.r_[alpha1, alphaN]
        
        factorial = N.array([reduce(lambda x,y:x*y, alphaN[0:n + 1]) for n in range(0, self.carbonNum)])
        sums = N.array([sum(factorial[0:n + 1]) for n in range(0, self.carbonNum)])
        sums = 1.0 + (1.0 + \
            pH2O / (self.K2 * self.K3 * self.K4 * pH2**2) + \
            1.0 / (self.K3 * self.K4 * pH2) + \
            1.0 / self.K4) * sums
        rCH4 = self.k5M * pH2 * alphaN[0] / sums[0]
        rParaffin = self.k5 * pH2 * factorial[1:] / sums[1:]
        rOlefin = self.k6 * (1.0 - betaN) * factorial[1:] / sums[1:]
        rCO2 = self.kv * (pCO * pH2O / pH2**0.5 - pCO2 * pH2**0.5 / self.kp) / \
            (1.0 + self.Kv * pCO * pH2O / pH2**0.5)
        # for C2+
        rCO = sum(N.arange(2.0, self.carbonNum + 1) * (rParaffin + rOlefin))
        rH2 = sum((2.0 * N.arange(2.0, self.carbonNum + 1) + 1.0) * rParaffin + \
            2.0 * N.arange(2.0, self.carbonNum + 1) * rOlefin)
        rH2O = rCO
        # modify for C1 and WGS
        rCO += rCH4 + rCO2
        rH2 += 3 * rCH4 - rCO2
        rH2O += rCH4 - rCO2
        return N.r_[-rCO, -rH2, rCH4, rH2O, rCO2, rParaffin, rOlefin]       
        
    
    def calc_heat(self, ra, sels):
        return 0.0
        

def get_ratelaw(type, params):
    """Factory function"""
    if type == 'L-H-2':
        return LH2RateLaw(params)
    elif type == 'DETAIL':
        return DetailRateLaw(params)
    else:
        raise ValueError('Unsupported rate law ID in the factory function')
         
         
        